Running the WBM code found on AWS:
####### Water Balance code ##########
source("WaterBalance/R/ET_functions.R")
source("WaterBalance/R/WB_functions.R")
############################################################# USER INPUTS #####################################################################
#KW comment: Historical was originally pulled from data we can't access. Therefore, I pulled this from AWS.
# (https://parkfutures.s3.us-west-2.amazonaws.com/maca-tprh-data/index.html)
Historical <- read.csv("katie_adds/centroid_data/PETE_historical.csv")
#Site characteristics
wb_sites <- read.csv("PETE_site_parameters.csv")
n<-nrow(wb_sites)
#Threshold temperature (deg C) for growing degree-days calculation
T.Base = 0
#Method for PET calculation
Method = "Oudin" #Oudin is default method for daily PRISM and MACA data (containing only Tmax, Tmin, and Date).
#Date format
DateFormat = "%m/%d/%Y"
############################################################ END USER INPUTS ###################################################################
############################################################ CREATE CLIMATE INPUTS #############################################################
#### Historical
# Convert pr.In to mm and F to C
Historical$ppt_mm <- (Historical$Precip..in.*25.4)
Historical$tmax_C <- 5/9*(Historical$Tmax..F. - 32)
Historical$tmin_C <- 5/9*(Historical$Tmin..F. - 32)
Historical$tmean_C <- 5/9*(Historical$Tavg..F. - 32)
######################################################### END CLIMATE INPUTS ####################################################################
######################################################### CALCULATE WB VARIABLES ################################################################
AllDailyWB<-list()
DailyWB<-Historical
for(i in 1:nrow(wb_sites)){
ID = wb_sites$SiteID[i] # KW: was previously "Site"
Lat = wb_sites$Lat[i]
Lon = wb_sites$Lon[i]
Elev = wb_sites$Elev[i] # KW: was previously "Elevation"
Aspect = wb_sites$Aspect[i]
Slope = wb_sites$Slope[i]
SWC.Max = wb_sites$SWC.Max[i]
Wind = wb_sites$Wind[i]
Snowpack.Init = wb_sites$Snowpack.Init[i]
Soil.Init = wb_sites$Soil.Init[i]
Shade.Coeff = wb_sites$Shade.Coeff[i]
#Calculate daily water balance variables
DailyWB$ID = ID
DailyWB$doy <- yday(DailyWB$Date)
DailyWB$daylength = get_daylength(DailyWB$Date, Lat)
DailyWB$jtemp = as.numeric(get_jtemp(Lon, Lat))
DailyWB$F = get_freeze(DailyWB$jtemp, DailyWB$tmean_C)
DailyWB$RAIN = get_rain(DailyWB$ppt_mm, DailyWB$F)
DailyWB$SNOW = get_snow(DailyWB$ppt_mm, DailyWB$F)
DailyWB$MELT = get_melt(DailyWB$tmean_C, DailyWB$jtemp, hock=4, DailyWB$SNOW, Snowpack.Init)
DailyWB$PACK = get_snowpack(DailyWB$jtemp, DailyWB$SNOW, DailyWB$MELT)
DailyWB$W = DailyWB$MELT + DailyWB$RAIN
if(Method == "Hamon"){
DailyWB$PET = ET_Hamon_daily(DailyWB)
} else {
if(Method == "Penman-Monteith"){
DailyWB$PET = ET_PenmanMonteith_daily(DailyWB)
} else {
if(Method == "Oudin"){
DailyWB$PET = get_OudinPET(DailyWB$doy, Lat, DailyWB$PACK, DailyWB$tmean_C, Slope, Aspect, Shade.Coeff)
} else {
print("Error - PET method not found")
}
}
}
DailyWB$PET = modify_PET(DailyWB$PET, Slope, Aspect, Lat, Shade.Coeff)
DailyWB$W_PET = DailyWB$W - DailyWB$PET
DailyWB$SOIL = get_soil(DailyWB$W, Soil.Init, DailyWB$PET, DailyWB$W_PET, SWC.Max)
DailyWB$DSOIL = diff(c(Soil.Init, DailyWB$SOIL))
DailyWB$AET = get_AET(DailyWB$W, DailyWB$PET, DailyWB$SOIL, Soil.Init)
DailyWB$W_ET_DSOIL = DailyWB$W - DailyWB$AET - DailyWB$DSOIL
DailyWB$D = DailyWB$PET - DailyWB$AET
DailyWB$GDD = get_GDD(DailyWB$tmean_C, T.Base)
AllDailyWB[[i]] = DailyWB
}
WBData<-do.call(rbind,AllDailyWB)
WBData_coded <- do.call(rbind,AllDailyWB) %>%
group_by(Date, GCM) %>%
dplyr::summarize(mean_runoff = mean(W_ET_DSOIL, na.rm = TRUE)/25.4,
mean_pet = mean(PET, na.rm = TRUE)/25.4,
mean_rain = mean(RAIN, na.rm = TRUE)/25.4) %>%
# To compare against our gridded product
filter(year(as_date(Date)) > 1979)
Pulling in then organizing the data associated with the gridded
product (pulling and averaging grids that intersect
wb_sites points) and centroid product on AWS.
# This is centroid data downloaded directly from AWS for PETE
# (https://parkfutures.s3.us-west-2.amazonaws.com/park-centroid-wb/index.html)
WBData_online_centroid <- read_csv("katie_adds/centroid_data/PETE_water_balance_historical.csv")%>%
filter(year(as_date(Date)) > 1979)
# This is the gridded data downloaded directly from AWS, the clipped to PETE's
# park boundary
WBData_online_grid <- rast("katie_adds/gridded_data/PETE_runoff_historic_gridmet_1980_2023.tif")
# next we clip the grid to only grids that intersect our wb_sites points in PETE:
NPS_points <- st_transform(sf::st_as_sf(wb_sites, coords = c("Lon", "Lat"), crs = 4326), st_crs(WBData_online_grid))
grid_runoff <- terra::extract(WBData_online_grid, NPS_points) %>%
pivot_longer(cols = -ID) %>%
dplyr::group_by(name) %>%
dplyr::summarize(mean_runoff = mean(value, na.rm = TRUE)/25.4/10) %>%
filter(name != "runoff_NA") %>%
mutate(date = as_date(sub(".*_", "", name))) %>%
filter(year(date) < 2023)
WBData_online_grid <- rast("katie_adds/gridded_data/PETE_PET_historic_gridmet_1980_2023.tif")
grid_pet <- terra::extract(WBData_online_grid, NPS_points) %>%
pivot_longer(cols = -ID) %>%
dplyr::group_by(name) %>%
dplyr::summarize(mean_pet = mean(value, na.rm = TRUE)/25.4/10) %>%
filter(name != "PET_NA") %>%
mutate(date = as_date(sub(".*_", "", name))) %>%
filter(year(date) < 2023)
WBData_online_grid <- rast("katie_adds/gridded_data/PETE_rain_historic_gridmet_1980_2023.tif")
grid_rain <- terra::extract(WBData_online_grid, NPS_points) %>%
pivot_longer(cols = -ID) %>%
dplyr::group_by(name) %>%
dplyr::summarize(mean_rain = mean(value, na.rm = TRUE)/25.4/10) %>%
filter(name != "rain_NA") %>%
mutate(date = as_date(sub(".*_", "", name))) %>%
filter(year(date) < 2023)
Here I am plotting each of the different WBM runoff products against each other for comparison.
year_colors <- c("#002EA3", "#C3F73A", "#E70870", "#256BF5", "#745CFB", "#FFCA3A", "#578010", "#56104E")
interpolated_colors <- colorRampPalette(year_colors)(length(unique(year(WBData_coded$Date))))
one <- ggplot() +
ggtitle("Runoff in") +
theme_bw() +
xlab("Amber's Code") +
ylab("Centroid Data on AWS") +
geom_point(aes(x = WBData_coded$mean_runoff,
y = WBData_online_centroid$`runoff in`,
color = factor(year(WBData_coded$Date)))) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
two <- ggplot() +
theme_bw() +
xlab("Amber's Code") +
ylab("Gridded Data on AWS") +
geom_point(aes(x = WBData_coded$mean_runoff,
y = grid_runoff$mean_runoff,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
three <- ggplot() +
theme_bw() +
xlab("Centroid Data on AWS") +
ylab("Gridded Data on AWS") +
geom_point(aes(x = WBData_online_centroid$`runoff in`,
y = grid_runoff$mean_runoff,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
ggpubr::ggarrange(one + theme(legend.position = "none"),
two + theme(legend.position = "none"),
three + theme(legend.position = "none"),
ncol = 1, nrow = 3,
common.legend = TRUE, legend = "right")
All three runoff products plotted together through time:
ggplot() +
geom_point(data = WBData_online_centroid, aes(x = as_date(Date), y = `runoff in`), color = "#E70870", alpha = 1) +
geom_point(data = grid_runoff, aes(x = date, y = mean_runoff), color = "#256BF5", alpha = 0.6) +
geom_line(data = WBData_coded, aes(x = as_date(Date), y = mean_runoff), color = "black", cex = 0.2) +
theme_bw() +
labs(x = "Date", y = "Runoff in") +
geom_text(data = WBData_online_centroid[1, ], aes(x = as.Date(Date), y = `runoff in`), label = "Centroid Data on AWS", color = "#E70870", vjust = -35, hjust = -1) +
geom_text(data = grid_runoff[1, ], aes(x = date, y = mean_runoff), label = "Gridded Data on AWS", color = "#256BF5", vjust = -33, hjust = -1) +
geom_text(data = WBData_coded[1, ], aes(x = as.Date(Date), y = mean_runoff), label = "Amber's Code", color = "black", vjust = -31, hjust = -1.8)
Here I am plotting each of the different WBM PET products against each other for comparison.
one <- ggplot() +
ggtitle("PET in") +
theme_bw() +
xlab("Amber's Code") +
ylab("Centroid Data on AWS") +
geom_point(aes(x = WBData_coded$mean_pet,
y = WBData_online_centroid$`PET in`,
color = factor(year(WBData_coded$Date)))) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
two <- ggplot() +
theme_bw() +
xlab("Amber's Code") +
ylab("Gridded Data on AWS") +
geom_point(aes(x = WBData_coded$mean_pet,
y = grid_pet$mean_pet,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
three <- ggplot() +
theme_bw() +
xlab("Centroid Data on AWS") +
ylab("Gridded Data on AWS") +
geom_point(aes(x = WBData_online_centroid$`PET in`,
y = grid_pet$mean_pet,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
ggpubr::ggarrange(one + theme(legend.position = "none"),
two + theme(legend.position = "none"),
three + theme(legend.position = "none"),
ncol = 1, nrow = 3,
common.legend = TRUE, legend = "right")
All three PET products plotted together through time:
ggplot() +
geom_point(data = WBData_online_centroid, aes(x = as_date(Date), y = `PET in`), color = "#E70870", alpha = 1) +
geom_point(data = grid_pet, aes(x = date, y = mean_pet), color = "#256BF5", alpha = 0.6) +
geom_line(data = WBData_coded, aes(x = as_date(Date), y = mean_pet), color = "black", cex = 0.2) +
theme_bw() +
labs(x = "Date", y = "PET in") +
geom_text(data = WBData_online_centroid[1, ], aes(x = as.Date(Date), y = `PET in`), label = "Centroid Data on AWS", color = "#E70870", vjust = -35, hjust = -1) +
geom_text(data = grid_pet[1, ], aes(x = date, y = mean_pet), label = "Gridded Data on AWS", color = "#256BF5", vjust = -33, hjust = -1) +
geom_text(data = WBData_coded[1, ], aes(x = as.Date(Date), y = mean_pet), label = "Amber's Code", color = "black", vjust = -31, hjust = -1.8)
Here I am plotting each of the different WBM rain products against each other for comparison.
one <- ggplot() +
ggtitle("Rain in") +
theme_bw() +
xlab("Amber's Code") +
ylab("Centroid Data on AWS") +
geom_point(aes(x = WBData_coded$mean_rain,
y = WBData_online_centroid$`rain in`,
color = factor(year(WBData_coded$Date)))) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
two <- ggplot() +
theme_bw() +
xlab("Amber's Code") +
ylab("Gridded Data on AWS") +
geom_point(aes(x = WBData_coded$mean_rain,
y = grid_rain$mean_rain,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
three <- ggplot() +
theme_bw() +
xlab("Centroid Data on AWS") +
ylab("Gridded Data on AWS") +
geom_point(aes(x = WBData_online_centroid$`rain in`,
y = grid_rain$mean_rain,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
ggpubr::ggarrange(one + theme(legend.position = "none"),
two + theme(legend.position = "none"),
three + theme(legend.position = "none"),
ncol = 1, nrow = 3,
common.legend = TRUE, legend = "right")
All three rain products plotted together through time:
ggplot() +
geom_line(data = WBData_coded, aes(x = as_date(Date), y = mean_rain), color = "black", cex = 0.2) +
geom_point(data = WBData_online_centroid, aes(x = as_date(Date), y = `rain in`), color = "#E70870", alpha = 1) +
geom_point(data = grid_rain, aes(x = date, y = mean_rain), color = "#256BF5", alpha = 0.6) +
#geom_line(data = WBData_coded, aes(x = as_date(Date), y = mean_rain), color = "black", cex = 0.2) +
theme_bw() +
labs(x = "Date", y = "Rain in") +
geom_text(data = WBData_online_centroid[1, ], aes(x = as.Date(Date), y = `PET in`), label = "Centroid Data on AWS", color = "#E70870", vjust = -35, hjust = -1) +
geom_text(data = grid_pet[1, ], aes(x = date, y = mean_pet), label = "Gridded Data on AWS", color = "#256BF5", vjust = -33, hjust = -1) +
geom_text(data = WBData_coded[1, ], aes(x = as.Date(Date), y = mean_pet), label = "Amber's Code", color = "black", vjust = -31, hjust = -1.8)